function map = glm_RFX_tMap(hfile, c, mapopts)
% GLM::RFX_tMap  - calculate a t contrast map
%
% FORMAT:       map = glm.RFX_tMap([c, mapopts])
%
% Input fields:
%
%       c           CxN contrast vector (default: full model and main eff)
%       mapopts     structure with optional fields
%        .interp    mesh-based interpolation (default: true)
%        .robust    use robust regression for estimation (default: false)
%        .srf       SRF object (required for multi-study MTC-based stats)
%
% Output fields:
%
%       map         MAP/VMP/SMP object with C maps

% Version:  v0.8a
% Build:    9111316
% Date:     Nov-13 2009, 4:54 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.ProjectTypeRFX ~= 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
numsubs = bc.NrOfSubjects;
numspred = bc.NrOfSubjectPredictors;
nval = numsubs - 1;
if nargin < 2 || ...
   ~isa(c, 'double') || ...
    isempty(c) || ...
    size(c, 1) > numspred || ...
    any(isinf(c(:)) | isnan(c(:)))
    c = [[eye(numspred - 1);zeros(1, numspred - 1)], [ones(numspred - 1, 1); 0]];
elseif numel(c) == size(c, 2)
    if numel(c) > numspred
        c((numspred + 1):end) = [];
    end
    c = real(c');
else
    c = real(c);
end
nummaps = size(c, 2);
if size(c, 1) < numspred
    c = [c; zeros(numspred - size(c, 1), size(c, 2))];
end
if nargin < 3 || ...
   ~isstruct(mapopts) || ...
    numel(mapopts) ~= 1
    mapopts = struct;
end
ipo = true;
if isfield(mapopts, 'interp') && ...
   ~isempty(mapopts.interp) && ...
   (isnumeric(mapopts.interp) || islogical(mapopts.interp))
    ipo = ipo && mapopts.interp(1);
end
if ~isfield(mapopts, 'robust') || ...
    numel(mapopts.robust) ~= 1 || ...
   ~islogical(mapopts.robust)
    mapopts.robust = false;
end

% generate map according to type
switch (bc.ProjectType)
    
    % currently unsupported types
    case {0}
        error( ...
            'BVQXfile:NotYetSupported', ...
            'FMR->MAP generation not yet supported.' ...
        );
    
    % VTCs
    case {1}
        szmap = size(bc.GLMData.RFXGlobalMap);
        map = BVQXfile('new:vmp');
        mapc = bvqxfile_getcont(map.L);
        mapc.NrOfMaps = size(c, 2);
        mapc.XStart = bc.XStart;
        mapc.XEnd = bc.XEnd;
        mapc.YStart = bc.YStart;
        mapc.YEnd = bc.YEnd;
        mapc.ZStart = bc.ZStart;
        mapc.ZEnd = bc.ZEnd;
        mapc.Resolution = bc.Resolution;
        mapc.Map.Type = 1;
        mapc.Map.LowerThreshold = -custom_tinv(0.005, nval);
        mapc.Map.UpperThreshold = -custom_tinv(0.0001, nval);
        mapc.Map.Name = 'Contrast:';
        mapc.Map.NrOfLags = [];
        mapc.Map.MinLag = [];
        mapc.Map.MaxLag = [];
        mapc.Map.CCOverlay = [];
        mapc.Map.ClusterSize = 25;
        mapc.Map.EnableClusterCheck = 0;
        mapc.Map.DF1 = nval;
        mapc.Map.DF2 = 0;
        mapc.Map.BonferroniValue = bc.NrOfVoxelsForBonfCorrection;
        mapc.Map.VMPData = single(zeros(szmap));
        mapc.Map(1:nummaps) = mapc.Map(1);
        
        % iterate over contrasts
        for cc = 1:size(c, 2)
            
            % get contrast
            conc = c(:, cc);
            coni = find(conc ~= 0);
            
            % initialize temp map
            tmpmp = zeros([szmap, numsubs]);
            for pc = coni(:)'
                for sc = 1:numsubs
                    tmpmp(:, :, :, sc) = tmpmp(:, :, :, sc) + conc(pc) * ...
                        bc.GLMData.Subject(sc).BetaMaps(:, :, :, pc);
                end
            end
            
            % generate 2nd-level single beta-t-maps
            if ~mapopts.robust
                tmap = sqrt(numsubs) * (mean(tmpmp, 4) ./ std(tmpmp, 0, 4));
            else
                [bmaps, wmaps] = fitrobustbisquare_img(ones(numsubs, 1), tmpmp);
                tmap = robustt(ones(numsubs, 1), tmpmp, bmaps, wmaps, 1);
            end
            tmap(isinf(tmap) | isnan(tmap)) = 0;
            
            % set name and map data
            mapc.Map(cc).Name = ...
                sprintf('Contrast:%s', sprintf(' %d', c(:, cc)'));
            mapc.Map(cc).VMPData = single(tmap);
        end
        
    % MTCs
    case {2}
        
        % SRF required
        numvert = bc.NrOfVertices;
        if ipo
            if ~isfield(mapopts, 'srf') || ...
                numel(mapopts.srf) ~= 1 || ...
               ~isBVQXfile(mapopts.srf, 'srf')
                error( ...
                    'BVQXfile:BadArgument', ...
                    'Missing or bad SRF reference in map options.' ...
                );
            end
            srfs = bvqxfile_getscont(mapopts.srf.L);
            srfc = srfs.C;
            if size(srfc.Neighbors, 1) ~= numvert
                error( ...
                    'BVQXfile:BadArgument', ...
                    'Number of vertices mismatch.' ...
                );
            end
            nei = srfc.Neighbors(:, 2);
            nnei = numel(nei);
            neil = zeros(nnei, 12);
            neis = zeros(nnei, 12);
            for nc = 1:nnei
                neis(nc) = numel(nei{nc});
                neil(nc, 1:neis(nc)) = nei{nc};
            end
            neii = (neil > 0);
        end
        
        map = BVQXfile('new:smp');
        mapc = bvqxfile_getcont(map.L);
        mapc.FileVersion = 4;
        mapc.NrOfVertices = numvert;
        mapc.NrOfMaps = size(c, 2);
        if ipo
            mapc.NameOfOriginalSRF = srfs.F;
        end
        mapc.Map.Type = 1;
        mapc.Map.NrOfLags = [];
        mapc.Map.MinLag = [];
        mapc.Map.MaxLag = [];
        mapc.Map.CCOverlay = [];
        mapc.Map.ClusterSize = 25;
        mapc.Map.EnableClusterCheck = 0;
        mapc.Map.DF1 = nval;
        mapc.Map.DF2 = 0;
        mapc.Map.BonferroniValue = bc.NrOfVoxelsForBonfCorrection;
        mapc.Map.Name = 'Contrast:';
        mapc.Map.SMPData = single(zeros(numvert, 1));
        mapc.Map(1:nummaps) = mapc.Map(1);
        
        % iterate over contrasts
        for cc = 1:size(c, 2)
            
            % get contrast
            conc = c(:, cc);
            coni = find(conc ~= 0);
            
            % initialize temp map
            tmpmp = zeros(numvert, numsubs);
            for pc = coni(:)'
                for sc = 1:numsubs
                    tmpmp(:, sc) = tmpmp(:, sc) + conc(pc) * ...
                        bc.GLMData.Subject(sc).BetaMaps(:, pc);
                end
            end
            
            % generate 2nd-level single beta-t-maps
            if ~mapopts.robust
                tmap = sqrt(numsubs) * (mean(tmpmp, 2) ./ std(tmpmp, 0, 2));
            else
                [bmaps, wmaps] = fitrobustbisquare_img(ones(numsubs, 1), tmpmp);
                tmap = robustt(ones(numsubs, 1), tmpmp, bmaps, wmaps, 1);
            end
            tmap(isinf(tmap) | isnan(tmap)) = 0;
            
            % interpolate ?
            if ipo
                rmap = zeros(numvert, 1);
                for nc = 1:size(neii, 2)
                    rmap(neii(:, nc)) = rmap(neii(:, nc)) + ...
                        tmap(neil(neii(:, nc), nc));
                end
                tmap = rmap ./ sum(neii, 2);
            end
            
            % set name and map data
            mapc.Map(cc).Name = ...
                sprintf('Contrast:%s', sprintf(' %d', c(:, cc)'));
            mapc.Map(cc).SMPData = tmap;
        end
        
     % unknown types
    otherwise
        error( ...
            'BVQXfile:UnknownOption', ...
            'Unknown GLM ProjectType: %d', ...
            bc.ProjectType ...
        );
end

% put back
bvqxfile_setcont(map.L, mapc);
